library("RColorBrewer")
library(Signac)
library(Seurat)
library(GenomicRanges)
library(future)
#library(SeuratWrappers)
library(harmony)
library(EnsDb.Hsapiens.v86)
library(stringr)
library(dplyr)
library(ggplot2)
library(patchwork)
library(kableExtra)
library(tidyverse)
set.seed(123)
# Paths
path_to_obj <- ("~/Documents/multiome_tonsil_Lucia/results/R_objects/13.tonsil_multiome_bcells_without_doublets_normalized.rds")
path_to_obj_all_data<- ("~/Documents/multiome_tonsil_Lucia/results/R_objects/12.tonsil_multiome_without_cluster7_doublets_normalized.rds")
path_to_markers<-("~/Documents/multiome_tonsil_Lucia/results/tables/tonsil_markers_bcell_075.csv")
tonsil_wnn_bcell <- readRDS(path_to_obj)
tonsil_wnn<-readRDS(path_to_obj_all_data)
tonsil_markers_075<-read_csv(path_to_markers)
## New names:
## * `` -> ...1
## Rows: 4688 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): ...1, gene
## dbl (6): p_val, avg_log2FC, pct.1, pct.2, p_val_adj, cluster
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
DimPlot(
tonsil_wnn_bcell,
group.by = "wsnn_res.0.075",
reduction = "wnn.umap",
pt.size = 0.1, label = T
)
tonsil_wnn_bcell$is_cluster4 <-
tonsil_wnn_bcell$wsnn_res.0.075 == "4"
tonsil_wnn_cluster4 <- subset(tonsil_wnn_bcell, subset = is_cluster4 == TRUE)
barcode_cluster4<-tonsil_wnn_cluster4@meta.data$lib_name_barcode
DimPlot(object = tonsil_wnn_bcell, cells.highlight = barcode_cluster4, reduction = "wnn.umap", cols.highlight = "red", cols = "gray", order = TRUE)+NoLegend()
DimPlot(object = tonsil_wnn_bcell, cells.highlight = barcode_cluster4, reduction = "umap.atac", cols.highlight = "red", cols = "gray", order = TRUE)+NoLegend()
DimPlot(object = tonsil_wnn_bcell, cells.highlight = barcode_cluster4, reduction = "umap.rna", cols.highlight = "red", cols = "gray", order = TRUE)+NoLegend()
joint.umap_age<- DimPlot(tonsil_wnn_bcell, label = FALSE, split.by = "age_group", pt.size = 0.1, reduction = "umap.rna") + plot_annotation(title = 'RNA UMAP ')+ ggtitle('age group')
joint.umap_hospital<- DimPlot(tonsil_wnn_bcell, label = FALSE, split.by = "hospital", pt.size = 0.1, reduction = "umap.rna") + plot_annotation(title = 'RNA UMAP ')+ ggtitle('hospital ')
joint.umap_age
joint.umap_hospital
joint.umap_age<- DimPlot(tonsil_wnn_bcell, label = FALSE, split.by = "age_group", pt.size = 0.1, reduction = "umap.atac") + plot_annotation(title = 'ATAC UMAP ')+ ggtitle('Age group')
joint.umap_hospital<- DimPlot(tonsil_wnn_bcell, label = FALSE, split.by = "hospital", pt.size = 0.1, reduction = "umap.atac") + plot_annotation(title = 'ATAC UMAP ')+ ggtitle('Hospital ')
joint.umap_age
joint.umap_hospital
DimPlot(object = tonsil_wnn, cells.highlight = barcode_cluster4, reduction = "wnn.umap", cols.highlight = "red", cols = "gray", order = TRUE)+NoLegend()
DimPlot(object = tonsil_wnn, cells.highlight = barcode_cluster4, reduction = "umap.atac", cols.highlight = "red", cols = "gray", order = TRUE)+NoLegend()
DimPlot(object = tonsil_wnn, cells.highlight = barcode_cluster4, reduction = "umap.rna", cols.highlight = "red", cols = "gray", order = TRUE)+NoLegend()
joint.umap_age<- DimPlot(tonsil_wnn, label = FALSE, split.by = "age_group", pt.size = 0.1, reduction = "umap.rna") + plot_annotation(title = 'RNA UMAP - All data ')+ ggtitle('age group')
joint.umap_hospital<- DimPlot(tonsil_wnn, label = FALSE, split.by = "hospital", pt.size = 0.1, reduction = "umap.rna") + plot_annotation(title = 'RNA UMAP ')+ ggtitle('hospital ')
joint.umap_age
joint.umap_hospital
joint.umap_age<- DimPlot(tonsil_wnn, label = FALSE, split.by = "age_group", pt.size = 0.1, reduction = "umap.atac") + plot_annotation(title = 'ATAC UMAP ')+ ggtitle('Age group')
joint.umap_hospital<- DimPlot(tonsil_wnn, label = FALSE, split.by = "hospital", pt.size = 0.1, reduction = "umap.atac") + plot_annotation(title = 'ATAC UMAP - All data ')+ ggtitle('Hospital ')
joint.umap_age
joint.umap_hospital
top10_tonsil_markers_075<-tonsil_markers_075 %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
top10mark_cluster4<-top10_tonsil_markers_075[["gene"]][41:50]
top7_tonsil_markers_075<-tonsil_markers_075 %>% group_by(cluster) %>% top_n(n = 7, wt = avg_log2FC)
top7mark_cluster4<-top7_tonsil_markers_075[["gene"]][29:35]
df_top10<-as.data.frame(top10_tonsil_markers_075)
DT::datatable(df_top10)
#install.packages("htmlwidgets", type = "binary")
#install.packages("DT", type = "binary")
df_top10<-as.data.frame(top10mark_cluster4)
df_top10
## top10mark_cluster4
## 1 FYB1
## 2 INPP4B
## 3 THEMIS
## 4 PRKCH
## 5 IL7R
## 6 LEF1
## 7 TC2N
## 8 CAMK4
## 9 BCL11B
## 10 ITK
#DoHeatmap(tonsil_wnn_bcell)
DoHeatmap(tonsil_wnn_bcell,top7_tonsil_markers_075$gene,assay = "RNA")
## Warning in DoHeatmap(tonsil_wnn_bcell, top7_tonsil_markers_075$gene, assay
## = "RNA"): The following features were omitted as they were not found in the
## scale.data slot for the RNA assay: FCER2, ZDHHC14
markers_kn<-c("BANK1", "FCER2","CD3D", "IL7R","MKI67", "TOP2A","MARCKSL1", "RGS13", "LMO2", "CCDC88A","GNLY", "NKG7","CD8A", "FCRL4", "FCRL5", "PLAC8", "SOX5","IGHG1", "IGHA1", "XBP1","LYZ", "S100A8", "CD3D")
DoHeatmap(tonsil_wnn_bcell,markers_kn,assay = "RNA")
## Warning in DoHeatmap(tonsil_wnn_bcell, markers_kn, assay = "RNA"): The following
## features were omitted as they were not found in the scale.data slot for the RNA
## assay: S100A8, PLAC8, CD8A, MARCKSL1, FCER2, BANK1
We are going to represent the top 7 markers
markers_gg_bcell <- function(x,umap){purrr::map(x, function(x) {
p <- FeaturePlot(
tonsil_wnn_bcell,
features = x,
reduction = umap,
pt.size = 0.1
)
p
})}
markers_gg_alldata <- function(x,umap){purrr::map(x, function(x) {
p <- FeaturePlot(
tonsil_wnn,
features = x,
reduction = umap,
pt.size = 0.1
)
p
})}
naive_mem_bcell<-c("BANK1", "FCER2")
cd4_tcell<-c("CD3D", "IL7R")
dz_gc_bcell<-c("MKI67", "TOP2A")
lz_gc_bcell<-c("MARCKSL1", "RGS13", "LMO2", "CCDC88A")
cytotoxic<-c("GNLY", "NKG7", "GZMK", "CD8A")
memory_bcell<-c( "FCRL4", "FCRL5", "PLAC8", "SOX5")
pc<-c("IGHG1", "IGHA1", "JCHAIN", "XBP1")
myeloid<-c("LYZ", "S100A8")
poor_q_doublets <-c("FDCSP", "CLU", "CXCL13", "CR2")
doublet_proliferative_tcell<-c("MT2A", "CD3D", "TRAC", "PCNA")
Unk<-c("PTPRCAP", "CD37", "CD74")
PDC<-c("PTCRA", "LILRA4", "IRF7")
DimPlot(
tonsil_wnn,
reduction = "wnn.umap",
pt.size = 0.1
) + ggtitle('WNN UMAP')
markers_gg_alldata(top7mark_cluster4,"wnn.umap")
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
DimPlot(
tonsil_wnn,
reduction = "umap.atac",
pt.size = 0.1
) + ggtitle('ATAC UMAP Harmony')
markers_gg_alldata(top7mark_cluster4,"umap.atac")
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
markers_gg_alldata(c("TCL1A","IGHD","FCER2","IGHM","IL4R"),"umap.atac")
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
markers_gg_alldata(naive_mem_bcell,"umap.atac")
## [[1]]
##
## [[2]]
markers_gg_alldata(memory_bcell,"umap.atac")
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
markers_gg_alldata(dz_gc_bcell,"umap.atac")
## [[1]]
##
## [[2]]
markers_gg_alldata(lz_gc_bcell,"umap.atac")
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
markers_gg_alldata(pc,"umap.atac")
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
markers_gg_alldata(cytotoxic,"umap.atac")
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
markers_gg_alldata(cd4_tcell,"umap.atac")
## [[1]]
##
## [[2]]
markers_gg_alldata(myeloid,"umap.atac")
## [[1]]
##
## [[2]]
DimPlot(
tonsil_wnn,
reduction = "umap.rna",
pt.size = 0.1
) + ggtitle('RNA Harmony')
markers_gg_alldata(top7mark_cluster4,"umap.rna")
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
markers_gg_alldata(c("TCL1A","IGHD","FCER2","IGHM","IL4R"),"umap.rna")
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
markers_gg_alldata(naive_mem_bcell,"umap.rna")
## [[1]]
##
## [[2]]
markers_gg_alldata(memory_bcell,"umap.rna")
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
markers_gg_alldata(dz_gc_bcell,"umap.rna")
## [[1]]
##
## [[2]]
markers_gg_alldata(lz_gc_bcell,"umap.rna")
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
markers_gg_alldata(pc,"umap.rna")
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
markers_gg_alldata(cytotoxic,"umap.rna")
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
markers_gg_alldata(cd4_tcell,"umap.rna")
## [[1]]
##
## [[2]]
markers_gg_alldata(myeloid,"umap.rna")
## [[1]]
##
## [[2]]
markers_gg_bcell(top7mark_cluster4,"wnn.umap")
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
DimPlot(
tonsil_wnn_bcell,
reduction = "umap.atac",
pt.size = 0.1
) + ggtitle('ATAC Harmony B cell')
markers_gg_bcell(top7mark_cluster4,"umap.atac")
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
DimPlot(
tonsil_wnn_bcell,
reduction = "umap.rna",
pt.size = 0.1
) + ggtitle('RNA Harmony B cell')
markers_gg_bcell(top7mark_cluster4,"umap.rna")
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] es_ES.UTF-8/es_ES.UTF-8/es_ES.UTF-8/C/es_ES.UTF-8/es_ES.UTF-8
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] forcats_0.5.1 purrr_0.3.4
## [3] readr_2.1.2 tidyr_1.2.0
## [5] tibble_3.1.6 tidyverse_1.3.1
## [7] kableExtra_1.3.4 patchwork_1.1.1
## [9] ggplot2_3.3.5 dplyr_1.0.8
## [11] stringr_1.4.0 EnsDb.Hsapiens.v86_2.99.0
## [13] ensembldb_2.18.2 AnnotationFilter_1.18.0
## [15] GenomicFeatures_1.46.1 AnnotationDbi_1.56.2
## [17] Biobase_2.54.0 harmony_0.1.0
## [19] Rcpp_1.0.8 future_1.23.0
## [21] GenomicRanges_1.46.1 GenomeInfoDb_1.30.0
## [23] IRanges_2.28.0 S4Vectors_0.32.3
## [25] BiocGenerics_0.40.0 SeuratObject_4.0.4
## [27] Seurat_4.1.0 Signac_1.5.0
## [29] RColorBrewer_1.1-2 BiocStyle_2.22.0
##
## loaded via a namespace (and not attached):
## [1] rappdirs_0.3.3 SnowballC_0.7.0
## [3] rtracklayer_1.54.0 scattermore_0.7
## [5] bit64_4.0.5 knitr_1.37
## [7] irlba_2.3.5 DelayedArray_0.20.0
## [9] data.table_1.14.2 rpart_4.1-15
## [11] KEGGREST_1.34.0 RCurl_1.98-1.5
## [13] generics_0.1.2 cowplot_1.1.1
## [15] RSQLite_2.2.9 RANN_2.6.1
## [17] bit_4.0.4 tzdb_0.2.0
## [19] spatstat.data_2.1-2 webshot_0.5.2
## [21] xml2_1.3.3 lubridate_1.8.0
## [23] httpuv_1.6.5 SummarizedExperiment_1.24.0
## [25] assertthat_0.2.1 xfun_0.29
## [27] hms_1.1.1 jquerylib_0.1.4
## [29] evaluate_0.14 promises_1.2.0.1
## [31] fansi_1.0.2 restfulr_0.0.13
## [33] progress_1.2.2 dbplyr_2.1.1
## [35] readxl_1.3.1 igraph_1.2.11
## [37] DBI_1.1.2 htmlwidgets_1.5.4
## [39] sparsesvd_0.2 spatstat.geom_2.3-1
## [41] ellipsis_0.3.2 crosstalk_1.2.0
## [43] backports_1.4.1 bookdown_0.24
## [45] biomaRt_2.50.1 deldir_1.0-6
## [47] MatrixGenerics_1.6.0 vctrs_0.3.8
## [49] ROCR_1.0-11 abind_1.4-5
## [51] cachem_1.0.6 withr_2.4.3
## [53] ggforce_0.3.3 vroom_1.5.7
## [55] sctransform_0.3.3 GenomicAlignments_1.30.0
## [57] prettyunits_1.1.1 goftest_1.2-3
## [59] svglite_2.0.0 cluster_2.1.2
## [61] lazyeval_0.2.2 crayon_1.5.0
## [63] pkgconfig_2.0.3 slam_0.1-49
## [65] labeling_0.4.2 tweenr_1.0.2
## [67] nlme_3.1-153 ProtGenerics_1.26.0
## [69] rlang_1.0.1 globals_0.14.0
## [71] lifecycle_1.0.1 miniUI_0.1.1.1
## [73] filelock_1.0.2 BiocFileCache_2.2.0
## [75] modelr_0.1.8 cellranger_1.1.0
## [77] polyclip_1.10-0 matrixStats_0.61.0
## [79] lmtest_0.9-39 Matrix_1.3-4
## [81] ggseqlogo_0.1 zoo_1.8-9
## [83] reprex_2.0.1 ggridges_0.5.3
## [85] png_0.1-7 viridisLite_0.4.0
## [87] rjson_0.2.20 bitops_1.0-7
## [89] KernSmooth_2.23-20 Biostrings_2.62.0
## [91] blob_1.2.2 parallelly_1.30.0
## [93] scales_1.1.1 memoise_2.0.1
## [95] magrittr_2.0.2 plyr_1.8.6
## [97] ica_1.0-2 zlibbioc_1.40.0
## [99] compiler_4.1.2 BiocIO_1.4.0
## [101] fitdistrplus_1.1-6 Rsamtools_2.10.0
## [103] cli_3.1.1 XVector_0.34.0
## [105] listenv_0.8.0 pbapply_1.5-0
## [107] MASS_7.3-54 mgcv_1.8-38
## [109] tidyselect_1.1.1 stringi_1.7.6
## [111] highr_0.9 yaml_2.2.2
## [113] ggrepel_0.9.1 grid_4.1.2
## [115] sass_0.4.0 fastmatch_1.1-3
## [117] tools_4.1.2 future.apply_1.8.1
## [119] parallel_4.1.2 rstudioapi_0.13
## [121] lsa_0.73.2 gridExtra_2.3
## [123] farver_2.1.0 Rtsne_0.15
## [125] digest_0.6.29 BiocManager_1.30.16
## [127] shiny_1.7.1 qlcMatrix_0.9.7
## [129] broom_0.7.12 later_1.3.0
## [131] RcppAnnoy_0.0.19 httr_1.4.2
## [133] colorspace_2.0-2 rvest_1.0.2
## [135] XML_3.99-0.8 fs_1.5.2
## [137] tensor_1.5 reticulate_1.24
## [139] splines_4.1.2 uwot_0.1.11
## [141] RcppRoll_0.3.0 spatstat.utils_2.3-0
## [143] plotly_4.10.0 systemfonts_1.0.3
## [145] xtable_1.8-4 jsonlite_1.7.3
## [147] R6_2.5.1 pillar_1.7.0
## [149] htmltools_0.5.2 mime_0.12
## [151] glue_1.6.1 fastmap_1.1.0
## [153] DT_0.20 BiocParallel_1.28.3
## [155] codetools_0.2-18 utf8_1.2.2
## [157] lattice_0.20-45 bslib_0.3.1
## [159] spatstat.sparse_2.1-0 curl_4.3.2
## [161] leiden_0.3.9 magick_2.7.3
## [163] survival_3.2-13 rmarkdown_2.11
## [165] docopt_0.7.1 munsell_0.5.0
## [167] GenomeInfoDbData_1.2.7 haven_2.4.3
## [169] reshape2_1.4.4 gtable_0.3.0
## [171] spatstat.core_2.3-2